library(psych)
library(ggplot2)
library(GGally)

EDA (Exploratory Data Analysis)

1. Прочитайте файл с данными

Загрузка данных

heart_data <- read.csv("/home/qgelena/code/scripts/R/stepik_projects/heart.csv")

Проверка структури дата фрейма

dim(heart_data)
## [1] 303  14
names(heart_data)
##  [1] "age"      "sex"      "cp"       "trestbps" "chol"     "fbs"     
##  [7] "restecg"  "thalach"  "exang"    "oldpeak"  "slope"    "ca"      
## [13] "thal"     "num"
head(heart_data)
##   age sex cp trestbps chol fbs restecg thalach exang oldpeak slope ca thal
## 1  63   1  1      145  233   1       2     150     0     2.3     3  0    6
## 2  67   1  4      160  286   0       2     108     1     1.5     2  3    3
## 3  67   1  4      120  229   0       2     129     1     2.6     2  2    7
## 4  37   1  3      130  250   0       0     187     0     3.5     3  0    3
## 5  41   0  2      130  204   0       2     172     0     1.4     1  0    3
## 6  56   1  2      120  236   0       0     178     0     0.8     1  0    3
##   num
## 1   0
## 2   2
## 3   1
## 4   0
## 5   0
## 6   0

2. Разберитесь, что значит каждый признак и в чем он измеряется

  1. age: age in years, непрерывный признак
  2. sex (1 = male; 0 = female), качественный признак
  3. cp: chest pain type, сердечно сосудистое заболевание, качественный – Value 1: typical angina – Value 2: atypical angina – Value 3: non-anginal pain – Value 4: asymptomatic
  4. trestbps: resting blood pressure (in mm Hg on admission to the hospital), непрерывный
  5. chol: serum cholestoral in mg/dl, сывороточный холестерин, непрерывный
  6. fbs: (fasting blood sugar > 120 mg/dl) (1 = true; 0 = false), сахар в крови натощак, качественный
  7. restecg: resting electrocardiographic results, качественный – Value 0: normal – Value 1: having ST-T wave abnormality (T wave inversions and/or ST elevation or depression of > 0.05 mV) – Value 2: showing probable or definite left ventricular hypertrophy by Estes’ criteria
  8. thalach: maximum heart rate achieved, непрерывный
  9. exang: exercise induced angina (1 = yes; 0 = no), качественный
  10. oldpeak = ST depression induced by exercise relative to rest, углубление/опускание ST во время нагрузки, сравнительно с отдыхом
  11. slope: the slope of the peak exercise ST segment, наклон ST сегмента во время пиковой нагрузки, качественный – Value 1: upsloping – Value 2: flat – Value 3: downsloping
  12. ca: number of major vessels (0-3) colored by flourosopy, дисретный
  13. thal (Thalium, a radioactive tracer injected during a stress test *): 3 = normal; 6 = fixed defect; 7 = reversable defect, качественный
  14. num: diagnosis of heart disease (angiographic disease status), качественный – Value 0: < 50% diameter narrowing – Value 1-4: > 50% diameter narrowing (in any major vessel: attributes 59 through 68 are vessels)

3. Проверьте пропущенные наблюдения

colSums(is.na(heart_data))
##      age      sex       cp trestbps     chol      fbs  restecg  thalach 
##        0        0        0        0        0        0        0        0 
##    exang  oldpeak    slope       ca     thal      num 
##        0        0        0        0        0        0

Работа с пропущенными значениями

Я решила удалить строки с пропущенными значениями, так как данных для анализа достаточно

heart_data$ca[heart_data$ca == "?"] <- NA
heart_data$thal[heart_data$thal == "?"] <- NA

colSums(is.na(heart_data))
##      age      sex       cp trestbps     chol      fbs  restecg  thalach 
##        0        0        0        0        0        0        0        0 
##    exang  oldpeak    slope       ca     thal      num 
##        0        0        0        4        2        0
heart_data <- na.omit(heart_data)

4. Определите тип каждого признака (количественные, порядковые, качественные).

Изменение типов переменных

Я решила изменить тип некоторых переменных и изменить их levels, так как в даном дата фрейме много векторов с качественными признаками

Качественные:

heart_data$sex <- factor(heart_data$sex)
levels(heart_data$sex) <- c("female", "male")

heart_data$cp <- factor(heart_data$cp)
levels(heart_data$cp) <- c("typical","atypical","non-anginal","asymptomatic")

heart_data$fbs <- factor(heart_data$fbs)
levels(heart_data$fbs) <- c("false", "true")

heart_data$restecg <- factor(heart_data$restecg)
levels(heart_data$restecg) <- c("normal","stt","hypertrophy")

heart_data$exang <- factor(heart_data$exang)
levels(heart_data$exang) <- c("no","yes")

heart_data$slope <- factor(heart_data$slope)
levels(heart_data$slope) <- c("upsloping","flat","downsloping")

heart_data$thal <- factor(heart_data$thal)
levels(heart_data$thal) <- c("normal","fixed","reversable")

heart_data$num <- ifelse(heart_data$num %in% c(0), 0, 1)

heart_data$num <- factor(heart_data$num)
levels(heart_data$num) <- c("Здоров", "Болен")

Дискретные: ca (не может количество сосудов быть дробным):

heart_data$ca <- factor(heart_data$ca) # не делаем конвертацию, так как это не обязательно;

Для количественных признаков определите, непрерывные они или дискретные

Лайфхак: если сомневаетесь, дискретный признак или нет, постройте его гистограмму с расзмером бина = 1. Если признак непрерывный, то одинаковые значения повторяются очень редко. Если дискретный, то вы увидите несколько высоких столбиков.

Я построила гистограммы для всех количественных признаков с размером бина 1, судя по выше сказаному, эти признаки непрерывные.

ggplot(heart_data, aes(x = heart_data$thalach))+geom_histogram(binwidth=1)

ggplot(heart_data, aes(x = heart_data$age))+geom_histogram(binwidth = 1)

ggplot(heart_data, aes(x = heart_data$trestbps))+geom_histogram(binwidth = 1)

ggplot(heart_data, aes(x = heart_data$chol))+geom_histogram(binwidth = 1)

ggplot(heart_data, aes(x = heart_data$oldpeak))+geom_histogram(binwidth = 1)

5. Смотрим на общие статистики: совпадают ли ожидаемые признаков значения с реальными?

Я использовала функцию describeBy, чтобы посмотреть на общие статистики, не учитывая данных с качественными признаками.

Можно заметить, что в выборке больше испытуемых, которые здоровы, 160 и 137, соответственно. В здоровых испытуемых можно увидеть крайние значения по переменной chol, сывороточный холестерин в отметке 564, что значительно больше, чем крайняя отметка для больных людей (409.0).

max кров. давление (thalach) в здоровых людей по среднему и медиане привышают значения в больных людей.

Также значения переменной oldpeak в здоровых людей ниже по медиане и среднему значению, чем в больных людей.

describe(heart_data)
##          vars   n   mean    sd median trimmed   mad min   max range  skew
## age         1 297  54.54  9.05   56.0   54.73  8.90  29  77.0  48.0 -0.22
## sex*        2 297   1.68  0.47    2.0    1.72  0.00   1   2.0   1.0 -0.75
## cp*         3 297   3.16  0.96    3.0    3.29  1.48   1   4.0   3.0 -0.84
## trestbps    4 297 131.69 17.76  130.0  130.48 14.83  94 200.0 106.0  0.69
## chol        5 297 247.35 52.00  243.0  244.70 47.44 126 564.0 438.0  1.11
## fbs*        6 297   1.14  0.35    1.0    1.06  0.00   1   2.0   1.0  2.01
## restecg*    7 297   2.00  0.99    2.0    2.00  1.48   1   3.0   2.0  0.01
## thalach     8 297 149.60 22.94  153.0  150.91 22.24  71 202.0 131.0 -0.53
## exang*      9 297   1.33  0.47    1.0    1.28  0.00   1   2.0   1.0  0.74
## oldpeak    10 297   1.06  1.17    0.8    0.88  1.19   0   6.2   6.2  1.23
## slope*     11 297   1.60  0.62    2.0    1.54  1.48   1   3.0   2.0  0.51
## ca*        12 297   1.68  0.94    1.0    1.51  0.00   1   4.0   3.0  1.17
## thal*      13 297   1.84  0.96    1.0    1.79  0.00   1   3.0   2.0  0.33
## num*       14 297   1.46  0.50    1.0    1.45  0.00   1   2.0   1.0  0.15
##          kurtosis   se
## age         -0.55 0.53
## sex*        -1.44 0.03
## cp*         -0.44 0.06
## trestbps     0.76 1.03
## chol         4.30 3.02
## fbs*         2.04 0.02
## restecg*    -1.99 0.06
## thalach     -0.09 1.33
## exang*      -1.46 0.03
## oldpeak      1.44 0.07
## slope*      -0.65 0.04
## ca*          0.19 0.05
## thal*       -1.83 0.06
## num*        -1.98 0.03
describeBy(heart_data[, -c(2, 3, 6, 7, 9, 11, 12, 13, 14)], heart_data$num, mat=T)
##           item group1 vars   n       mean         sd median     trimmed
## age1         1 Здоров    1 160  52.643750  9.5511512   52.0  52.4531250
## age2         2  Болен    1 137  56.759124  7.8996701   58.0  57.2432432
## trestbps1    3 Здоров    2 160 129.175000 16.3739900  130.0 128.5390625
## trestbps2    4  Болен    2 137 134.635036 18.8967302  130.0 133.1171171
## chol1        5 Здоров    3 160 243.493750 53.7575499  235.5 239.2421875
## chol2        6  Болен    3 137 251.854015 49.6799374  253.0 251.0000000
## thalach1     7 Здоров    4 160 158.581250 19.0433045  161.0 160.0468750
## thalach2     8  Болен    4 137 139.109489 22.7106735  142.0 139.7837838
## oldpeak1     9 Здоров    5 160   0.598750  0.7871601    0.2   0.4617188
## oldpeak2    10  Болен    5 137   1.589051  1.3050061    1.4   1.4747748
##                mad min   max range       skew   kurtosis         se
## age1      10.37820  29  76.0  47.0  0.1191477 -0.6652838 0.75508480
## age2       5.93040  35  77.0  42.0 -0.5736352  0.1352938 0.67491436
## trestbps1 14.82600  94 180.0  86.0  0.4346036  0.2568030 1.29447757
## trestbps2 14.82600 100 200.0 100.0  0.8014577  0.6514438 1.61445661
## chol1     44.47800 126 564.0 438.0  1.7032585  7.2146181 4.24990748
## chol2     50.40840 131 409.0 278.0  0.2771550  0.2207174 4.24444349
## thalach1  16.30860  96 202.0 106.0 -0.6789875  0.4109986 1.50550541
## thalach2  25.20420  71 195.0 124.0 -0.2857344 -0.3011150 1.94030378
## oldpeak1   0.29652   0   4.2   4.2  1.5772389  2.8336782 0.06223047
## oldpeak2   1.48260   0   6.2   6.2  0.7130820  0.2774105 0.11149420

6. Строим scatterplot для всех пар, долго его разглядываем с точки зрения outliers, неоднородностей, вида распределений, вида зависимостей (линейные/нелинейные) и прочие. Запишите свои наблюдения.

Строим scatterplot для всех пар

pairs(heart_data[, -c(2, 3, 6, 7, 9, 11, 12, 13, 14)], col=heart_data$num)

ggpairs(heart_data[, -c(2, 3, 6, 7, 9, 11, 12, 13, 14)])

#Но в таком виде совершенно ничего не видно :( 
#Советую узнать, как изменить размер графика и сделать его больше
ggpairs(heart_data, lower = list(combo = wrap("facethist", binwidth = 1, mapping=aes(colour=num, alpha=0.4))))

7. Если есть outliers, то попробуйте объяснить их причину (ошибка в данных, особые индивиды) и удалить их

Было замечено несколько аутлаеров, некоторые из них при этом находятся в категории “здоровые”. Например, женщина с chol больше 500 и thalach около 160. Я бы уточнила у авторов исследования, не являются ли подобные значения ошибками измерений

pairs(heart_data[ ,c(1,4,5,8,10)], col = heart_data$num)
par (xpd=TRUE)
legend("bottomright", fill=unique(heart_data$num), legend=c(levels(heart_data$num)))

8. Если есть неоднородности (например, видны два облака точек), то постарайтесь объяснить причину (найти категоризующую переменную, объясняющую эту неоднородность). Для этого нужно раскрасить scatterplot для всех пар в цвета по всем категориальным переменным и найти ту, такую раскраску, в которой ваши облака одноцветные.

Хорошо заметна разница между группами здровых и больных по показателям oldpeak, thalach. У здоровых выше показатель thalach, у больных выше oldpeak

pairs(heart_data[ ,c(1,4,5,8,10)], col = heart_data$slope)

pairs(heart_data[ ,c(1,4,5,8,10)], col = heart_data$cp)

Можно заметить, что аутлаеры по переменной oldpeak лежат в одной (зелёной) категории slope. Также они лежат в одной (синей) категории cp

Разделения данных на чёткие облака отмечено не было

#“аутлайнеры по переменной oldpeak лежат в одной (зелёной) категории slope” а какому значению соответсвует зеленый slope?

pairs(heart_data[ ,c(1,4,5,8,10)], col = heart_data$slope)
par (xpd=TRUE)
legend("bottomright", fill=unique(heart_data$slope), legend=c(levels(heart_data$slope)))

9. Графики с помощью ggplot2. Минимум 5 разнообразных графиков. Обязательно постройте эти графики и опишите словами, что интересного вы в них нашли.

гистограммы (не забудьте поэкспериментировать с размером бинов и объяснить, как изменился вид гистограммы) ящики с усами и с точечками вайлин плоты с точечками попробуйте facet_grid разбивайте графики по категориальным переменным и раскрашивайте их в разные цвета (попробуйте разные цветовые схемы)

ggpairs(heart_data[ ,c(1,4,5,8,10)], aes(col=heart_data$num, alpha = 0.4))

У здоровых выше показатель thalach, у больных выше oldpeak

ggpairs(heart_data[ ,c(1,4,5,8,10)], aes(col=heart_data$slope, alpha = 0.4))

“зелёный” slope - это “2” Value 2: flat" - плоский угол

ggpairs(heart_data, aes(col=heart_data$num, alpha = 0.4))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(heart_data, aes(x=num, y=oldpeak)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=4) +
  geom_jitter(shape=16, position=position_jitter(0.2))

На этом графике мы видим, что по показателю oldpeak у больных, например, существенно больше медиана, а также “намного больше сам ящик”, то есть выше среднеквдратическое отклонение (Надеюсь, я не путаюсь в том, что ширина box plot позволеят судить о среднеквадратическом отклонении, раз для его построения мы откладываем 1,5 квартиля, размер которых как раз зависит от sd)

#{r} #ggplot(heart_data, aes(x=age, y=oldpeak, col=thalach)) + geom_point() + facet_grid(vars(num), vars(sex)) #

На этом графике видно, что больных женщин в испытании было существенно меньше. Если игнорировать outliers, то можно сказать, что характер распределения признака среди больных и здоровых похож у обоих полов

ggplot(heart_data, aes(x=num, y=oldpeak, col=sex)) + 
  geom_violin() +
  geom_jitter(shape=16, position=position_jitter(0.2))

Различия в распределении признака идут по параметру здоровый/больной, а не полу, что различаются медианы и отклонение

ggplot(heart_data, aes(thalach, fill = sex, alpha = 0.5)) +
  geom_histogram(binwidth = 5)

На этой гистограмме мы видим, что у женщин показатель thalach в среднем выше

ggplot(heart_data, aes(thalach, fill = sex, alpha = 0.5)) +
  geom_histogram(binwidth = 50)

Давайте поэкспериментируем с bin`ом. Можно заметить, что при неправильно подобранном шаге информативность гистограммы существенно снижается

Графики

#я бы еще добавила легенду к цветам, чтобы было еще понятнее
ggplot(data=heart_data[!is.na(heart_data$thalach),], aes(num, thalach)) + 
  geom_violin(aes(fill = num), trim = FALSE, alpha = 0.3) + 
  geom_boxplot(aes(fill = num), width = 0.2, outlier.colour = NA) + 
  theme(legend.position = "NA") + geom_jitter(shape=16, position=position_jitter(0.2))

heart_data$new.var <- heart_data$age > 50
ggplot(heart_data, aes(x=chol, y=slope)) +
geom_point(alpha=0.5) +
facet_grid(new.var ~sex) +
scale_color_gradient(low = "blue", high = "red") +
theme_bw()

ggplot(heart_data, aes(x = chol)) + geom_histogram(fill='blue', col="yellow", binwidth = 7)

hist(heart_data$chol)

# Так никогда не делать (!):
# ggplot(heart_data, aes(x=oldpeak, fill =slope)) + geom_dotplot()

ggplot(heart_data, aes(x=num, y=chol)) + geom_boxplot()

Review

  1. В чанках, которые кидают ворнинги используй waning =F и message = F, чтобы не захламлять html. Например, вот тут стоило отключить ворнинги {r} library(psych) library(ggplot2) library(GGally)

  2. Про количественные переменные нужно было понять дискретные они или непрерывные.

  3. Вот эта операция совсем некорректная heart_data_dis <- heart_data[heart_data$num == "0" | heart_data$num == "1", ] 0 кодирует здоровых людей, а 1,2,3,4 больных. Тебе нужно было склеить всех больных людей, а ты вместо этого выкинула тех, у кого heart_data\(num == "2" или heart_data\)num == “3” или heart_data$num == “4”, то есть выкинула большой кусок датасета.

  4. Там, где рисуешь парные скаттерплоты нужно описывать, что ты на них видишь – зависимости, неоднородности в данных, аутлаеры, …

  5. Здорово, что научилась строить такой график ggpairs(heart_data_dis, lower = list(combo = wrap(“facethist”, binwidth = 1, mapping=aes(colour=num, alpha=0.4)))) Но в таком виде совершенно ничего не видно :( Советую узнать, как изменить размер графика и сделать его больше

  6. Вот этот график классный! ggplot(data=heart_data_dis[!is.na(heart_data_dis$thalach),], aes(num, thalach)) + geom_violin(aes(fill = num), trim = FALSE, alpha = 0.3) + geom_boxplot(aes(fill = num), width = 0.2, outlier.colour = NA) + theme(legend.position = “NA”) + geom_jitter(shape=16, position=position_jitter(0.2)) я бы еще добавила легенду к цветам, чтобы было еще понятнее